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NONPARAMETRIC TESTS OF STRUCTURE FOR HIGH 
ANGULAR RESOLUTION DIFFUSION IMAGING IN Q-SPACE 1 

By Sofia C. Olhede and Brandon Whitcher 

University College London and Glaxo SmithKline 

High angular resolution diffusion imaging data is the observed 
characteristic function for the local diffusion of water molecules in 
tissue. This data is used to infer structural information in brain 
imaging. Nonparametric scalar measures are proposed to summarize 
such data, and to locally characterize spatial features of the diffu- 
sion probability density function (PDF), relying on the geometry of 
the characteristic function. Summary statistics are defined so that 
their distributions are, to first-order, both independent of nuisance 
parameters and also analytically tractable. The dominant direction 
of the diffusion at a spatial location (voxel) is determined, and a new 
set of axes are introduced in Fourier space. Variation quantified in 
these axes determines the local spatial properties of the diffusion 
density. Nonparametric hypothesis tests for determining whether the 
diffusion is unimodal, isotropic or multi-modal are proposed. More 
subtle characteristics of white-matter microstructure, such as the de- 
gree of anisotropy of the PDF and symmetry compared with a variety 
of asymmetric PDF alternatives, may be ascertained directly in the 
Fourier domain without parametric assumptions on the form of the 
diffusion PDF. We simulate a set of diffusion processes and charac- 
terize their local properties using the newly introduced summaries. 
We show how complex white- matter structures across multiple voxels 
exhibit clear ellipsoidal and asymmetric structure in simulation, and 
assess the performance of the statistics in clinically-acquired magnetic 
resonance imaging data. 

1. Introduction. Many applications in brain imaging are based on cal- 
culating local statistics that are later combined to infer global properties 
of spatial links or functional connections. In this paper we focus on the 
local analysis of high angular resolution diffusion imaging (HARDI) data, 
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a special type of magnetic resonance imaging (MRI). HARDI observations 
correspond to the local (in a single voxel 2 ) measurement of the local molec- 
ular diffusion of water at a number of different orientations over a spherical 
shell of fixed radius [Callaghan (1993)]. Measurements from an MRI scan- 
ner are taken directly in the Fourier domain and translated into the spatial 
domain via the inverse Fourier transform. 

A HARDI acquisition scheme permits the characterization of directional 
spatial properties of the diffusion probability density function (PDF). The 
local structure of white-matter brain tissue may be inferred from such mea- 
surements [Basser, Mattiello and Bihan (1994); Basser (2002)]. Once local 
statistics have been formed, it is of interest to combine information across 
voxels (spatial locations), for example, to connect local directions of esti- 
mated diffusion PDFs to recognize major nerve fiber tracts, to infer local 
fiber structure from the estimated diffusions [Mori and van Zijl (2002)], 
and/or to use other locally-defined statistical summaries in inferential pro- 
cedures [Jensen et al. (2005)]. 

Different orientational sampling designs can be used at each voxel and, if 
a simple parametric model is used for the PDF, then rather sparse sampling 
will be sufficient to recover the parameters of the model. Traditional analysis 
of HARDI measurements is based on modeling the diffusion PDF paramet- 
rically as a (zero- mean) Gaussian, and estimating a diffusion tensor (the 
covariance matrix of the Gaussian PDF), a procedure which corresponds 
to diffusion tensor imaging (DTI). Such methods have drawbacks, namely, 
of not describing more complex white-matter structures well, and their us- 
age trades a small variance for potentially large bias. While the diffusion 
tensor model has both theoretical justification — and has been extremely 
popular — it prohibits one from describing more complicated white-matter 
microstructure, such as crossing, kissing and forking fibers [Mori and van 
Zijl (2002)]. 

It is believed that intravoxel orientational heterogeneity affects as many 
as one third of all imaged white-matter voxels [Behrens et al. (2007)], and so 
addressing such structure is important. With more time- intensive sampling 
schemes (such as HARDI [Tuch et al. (2002)] or diffusion spectrum imag- 
ing), the possibility of more complicated estimators may be used, for exam- 
ple, multi-tensor modeling [Alexander (2005)], nonparametric alternatives 
such as persistent angular structure MRI [Jansons and Alexander (2003)], 
Q-ball imaging [Tuch (2004)], the diffusion orientation transform [Ozarslan 
et al. (2006)] and spherical deconvolution [Tournier et al. (2004)]. While 
using a nonparametric approach removes bias, usage of such nonparamet- 
ric methods is challenging because the diffusion process is measured in the 



2 A voxel is a three-dimensional "volume element" of data, just as a pixel is a two- 
dimensional "area element" of data. 
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Fourier domain (q-space 3 ), and the characteristic function has been consid- 
erably undersampled to accommodate realistic scanning times in practice. 
This challenges the stable inversion of information, the local characteristic 
function, to local spatial structure. 

This paper develops a statistical framework, using nonparametric meth- 
ods, for characterizing HARDI data directly in g-space [Tuch et al. (2002)] 
without local inversion. This avoids calculating nonlinear transformations of 
the data, whose usage usually leads to intractability of the distributions of 
statistical summaries. The approximate distributions of the proposed esti- 
mators in this paper are derived and are defined so that, to first order, they 
are free of any nuisance parameters. The proposed statistics are a first step 
toward the automated detection of subtle characteristics of white-matter 
microstructure, that is, scalene diffusions (Figure 1) or asymmetry in de- 
cay in a fixed axis. Both properties, scalene diffusion and asymmetry, have 
been found in a forking fiber structure (Figure 1), and may be important 
summaries to feed into fiber-tracking algorithms [Mori and van Zijl (2002)]. 
The derived methods also serve as a warning when interpreting multi-tensor 
models in clinically-feasible acquisition schemes, as similar characteristics 
can be obtained from more complex single peaked structures. 

Global features like bi- or multi-modality of the diffusion PDF are de- 
scribed reasonably well by many methods over a range of signal-to-noise 
ratios (SNRs), with the small caveat that the various implicit assumptions 
inherent to any of the given methods must be satisfied. Parametric models 
introduce bias when they are not appropriate, whereas using a nonpara- 
metric method increases the variance in the estimation. Using a moderate 
number of directions in the HARDI sampling scheme restricts the possi- 
bility of determining smaller scale structure of the diffusion PDF. Strong 
parametric assumptions increase the power of any proposed statistic to de- 
tect multiple diffusion directions, with the consequence that any deviation 
from the prescribed structure in the parametric model may be used to reject 
null hypotheses such as unimodality. 

In the method proposed here to determine the properties of the diffu- 
sion PDF, prolate diffusion PDFs are separated from isotropic (or spher- 
ical) PDFs using a test based on a comparison of relative magnitudes in 
(/-space; see Figure 1 for illustrations of prolate and spherical diffusion mod- 
els. Subsequently, multi-modal distributions are then differentiated from the 
isotropic and unidirectional. The unidirectional diffusion is associated with 



3 Q-space is the Fourier domain representation of the local diffusion and is the space 
where measurements are made in MRI. The global image Fourier representation is usually 
inverted to a spatial representation, but the local Fourier transform is not inverted as part 
of the acquisition, leaving the spatial domain observations associated with a measurement 
of local diffusion in a Fourier domain orientation. 
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Fig. 1. Simplified diagrams for typical Gaussian diffusion models (first column) and 
fiber configurations in a voxel of white matter in the brain (second column). Spherical 
diffusion is found when no fibers are present in a voxel of brain tissue (e.g., cerebral-spinal 
fluid) and all eigenvalues are equal (Ai = A2 = A3). Prolate diffusion is when a single fiber 
bundle is present in the voxel (Ai 2> A2 = A3). Scalene diffusion is when two fiber bundles 
of similar mass cross in perpendicular directions (Ai ~ A2 3> A3) . The concept of "crossing 
fibers" involves two fiber bundles that do not necessarily intersect at right angles in the 
same voxel. The concept of "kissing fibers" involves two fiber bundles that occupy the same 
voxel, but do not intersect. The concept of "forking fibers" involves a single fiber going in 
the voxel and two fiber bundles leaving the voxel. A "fanning fiber" (not shown) is similar 
to a forking fiber, but instead of a single direction the fiber produces multiple diverging 
fibers on one side of the voxel. 



a great circle in (/-space [Tuch (2004)], and we call this the dominant great 
circle. The strongest direction defines an important spatial summary of the 
diffusion PDF, and specifies the major axis of the diffusion in g-space (Fig- 
ure 2g). The perpendicular to the major direction in space defines a set of 
points lying on a great circle in g-space, which exactly corresponds to the 
dominant great circle. 

If a given voxel has been diagnosed as unidirectional (or if there is a dom- 
inant great circle in q-space), then we seek to characterize its main unidi- 
rectional structure in more detail. A scalar measure resembling the popular 
fractional anisotropy 4 is defined as the anisotropy statistic, by comparing the 



4 The fractional anisotropy (FA) is a measure of uniformity of the eigenvalues of a Gaus- 
sian covariance matrix [Basser and Pierpaoli (1996)]. 
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magnitude of the g-space diffusion on the dominant great circle with its two 
perpendicular point (s). This measure determines the degree of anisotropy 
of the diffusion PDF. Further investigation of unidirectional voxels causes 
us to focus on quantifying the uniformity of decay in the minor axes of the 
diffusion PDF, or the perpendicular to the dominant great circle, to describe 
further detailed structure of the characteristic function. 

Ellipsoidal diffusions are an important class of diffusions and the sca- 
lene structure of the diffusion PDF is particularly important when com- 
bining voxel- wise information [Seunarine et al. (2007)]. The aforementioned 
work showed that the scalene structure of the peak is related to the peak 
anisotropy in space and important for treating bending and fanning fibers 
(Figure 1). For diffusions with ellipsoidal decay, their minor axes are well 
defined by this (scalene) decay structure, while for nonellipsoidal diffusions 
the minor axes correspond to a set of axes in the plane of the dominant 
great circle, parameterizing locations on the dominant great circle. We ex- 
amine the scalene structure of the diffusion PDF, which is quantified by 
the difference in decay in the two spatial minor axes, defined as such also 
for nonellipsoid diffusions. This corresponds to examining the variability of 
the diffusion on the great circle perpendicular to the vector associated with 
the major direction of the diffusion. For a Gaussian diffusion model this is 
given by the two minor eigenvalues of the eigen-decomposition of the diffu- 
sion tensor. A statistical test for uniformity on the great circle is developed 
that can be related to the spatial decay of the diffusion PDF in the minor 
axes. Another feature of interest in the PDF is asymmetry in the decay in 
a fixed direction perpendicular to the dominant great circle. This heuris- 
tic may be visualized in space as a diffusion PDF that appears ellipsoidal 
but the peak is in one of the foci rather than the center of the ellipse. We 
introduce a test statistic for asymmetry based on this understanding. To 
motivate our interest in asymmetry and ellipsoidality, we simulate forking 
and crossing structures, and show how both asymmetry and ellipsoidality 
follow as precursors to forking structure, and such information could be used 
to improve the tracking of fibers. 

The methodology presented here improves our understanding of the dif- 
fusion PDF by not relying on parametric assumptions when analyzing the 
measurements, yet still relating (/-space structure directly to spatial proper- 
ties. Nonparametric statistical summaries are defined directly in q-space to 
increase the power of the proposed hypothesis tests and theoretical critical 
values for the statistics are provided. Understanding the inherent limita- 
tions of HARDI measurements can be obtained directly from our discussion 
of simulated diffusions, thus increasing the understanding of parametric as- 
sumptions that are necessary to derive more complicated structures from 
the diffusion PDF. 
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2. Statistical models for HARDI data. 

2.1. Observational model. We denote the sampling of the observations 
by the set Qo = {«li}f=i- At each qj = (qn,qi2,Qi3) on the unit sphere ||q|| = 
(of + <?2 + ^i) 1 ^ 2 = 1 we obtain an observed measurement -A(qi) > 0, corre- 
sponding to the magnitude of a complex- valued observation (proportional to 
the noisy characteristic function of local diffusion 5 ). Furthermore, we take no 
observations at q = 0, denoted by Ak(0) for k = 1, . .. , no- We distinguish 
here between the measured apparent diffusion at qj, namely, A{c\j), and the 
theoretical diffusion value, A(c[i). Note that the expected value of A(c[i) is 
not equivalent to A{<\i), for two reasons. First because the observations are 
magnitudes, with the noise contributing in the expectation, and second we 
need to re-normalize the observed diffusion to have unit volume, as noted by 
Alexander (2005). As the PDF is a density, it has to satisfy the normalization 
of 

(2.1) JJJa(x)d 3 x = l .4.(0) = 1, 

where a(x) is the diffusion probability density function (PDF), or the in- 
verse Fourier Transform of »4(q). We apply a biased estimator of a sim- 
ple average to estimate the inverse of the normalizing constant by A(0) = 
n o 1 ^2k=i Ak(0)- We re-normalize the observed diffusion such that -A(qi) = 
A((\i)/A(0). The diffusion value A(q^) has (approximately) a Rician distri- 
bution with parameters «4(qj) and a 2 [Gudbjartsson and Patz (1995)]. As 
the SNR will be large at q = 0, the noise floor of the Rician distribution will 
have limited impact in the estimation of the normalization constant. While 
the diffusion PDF a(x) is not Gaussian, the Rician distribution under reason- 
able SNR is well approximated by the Gaussian, and sums of Rician variables 
will be very similar to a Gaussian. In subsequent sections we shall calculate 
statistical estimators from normalized measurements {^(qj)}™ and look at 
maxima of these statistics, which may be represented (approximately) by the 
maxima of suitably-scaled Gaussian random variables. If we are in the regime 
of low SNR, then these test statistics will be approximated by a mixture of 
Gaussian and Chi random variables whose tail-behavior is not substantially 
heavier than Gaussian random variables, but whose mean is not consistent 
with our results. An assumption for the method to work is therefore a rea- 
sonable level of the SNR, as is further discussed in Section 5. 

The normalized diffusion measurements A(c[i) should exhibit symme- 
try as the diffusion PDF is real- valued, symmetric and indeed positive, 



5 Note that this is different from the empirical characteristic function. 
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that is, A(-q.i) = A(qi) [Wedeen et al. (2005)]. To fully exploit the Her- 
mitian symmetry, we shall reflect the observations to the augmented set 
Q = {q:q G Qo} U {q:-q G Qo}, and set A{— q 4 ) = A(q,) [Jansons and 
Alexander (2003)]. 

We assume that a nonparametric estimator of the diffusion in g-space is 
constructed. For our purposes we have chosen to use a variable-bandwidth 
estimator [Olhede and Whitcher (2008a, 2008b)], but the methodology out- 
lined here is applicable to other linear estimators (e.g., radial basis functions 
and /or spherical harmonics) with some straightforward alteration of the sta- 
tistical properties (specifically, second-order structure) of the estimators. 

2.2. Great circles in q-space. Spatial properties of the diffusion PDF may 
be described directly in g-space. The advantage of such an operation is that 
we avoid the need to invert the PDF to the spatial domain for analysis, al- 
lowing us to employ a broad range of modeling approaches. A basic building 
block of our analysis is an ellipsoid density. We refer to a density cie(x) as 
an ellipsoid density if its FT takes the form 



(2.2) A E (^A,T) = B 

where Xj > for j = 1,2,3, {vj} constitutes a basis for ]R 3 and B(-) is 
a monotonically decreasing function. For example, it is common to use the 
Gaussian characteristic function B(q) = e~ 2 ^ nq ^ . We collect the eigenvalues 
in the matrix A = diag(Ai, A2, A3), and define 

vu V12 V13 
V21 v 2 2 v 2 3 

V31 V 32 V 33 ,_ 

to model the axis of any orientational structure. Ellipsoid densities are nat- 
ural building blocks, just like the special case of the DTI model, but do not 
(for example) include multi-modal densities. If the (/-space density takes this 
form, then the spatial PDF is given by inverting the FT 

(2.4) a E (x; A, T) = /// .4 E (q; A, r)e i2 ^ d 3 q 

[Callaghan (1993)]. We note for x G M 3 , with x = ||x|| and q = ||q||, that 
oe(x;A,T) takes the form 

(2.5) «e(x; A,T) = |A| 1 /2 6 (|| A - 1 /2 Tx ||) ) 



E 



q| 



(2.3) 



-T 
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where 

/oo roo roo 
/ / B(qy 27TqT *d 3 q 
-oo J —oo J —oo 

(2.7) -ikr^)***** 

2 r°° 

(2.8) =- / B(q)sm(2Trxq)qdq, 

x Jo 

which follows from Gradshteyn and Ryzhik (2000), page 1112. The meaning 
of "ellipsoid density" becomes clear from this expression, since whenever 
||A _1 / 2 Tx|| = R, where R > is a constant, the function «e( x ) takes the 
same value in space. As long as all the eigenvalues are positive, cie(-) will 
map out ellipsoidal contours of equal function value in space. The Gaussian 
DTI model fits into this class of densities with b(x) = (27r)~ 3 ^ 2 e~ x I 2 as 
well as, for example, the Matern family with the spatial variable exchanged 
with the spatial-frequency variable [Matern (I960)]. The model proposed by 
Kaden, Knosche and Anwander (2007) is also related to such densities. 

Figure 2 provides examples of diffusion processes displayed in both the 
spatial and frequency ((/-space) domains. The spatial domain corresponds to 
the diffusion PDF, whereas its Fourier transform corresponds to the (/-space 
representation. Common processes, such as prolate and scalene diffusion, are 
given as well as more exotic examples, such as a mixture of prolate diffusion 
processes and a process that cannot be represented using a Gaussian diffu- 
sion model. The values of T specify the orientation of the diffusion PDF, 
while A gives its qualitative appearance when coupled with B(-). Looking 
directly at Figure 2, it may be difficult for one to appreciate the local struc- 
ture near the peak, which motivates us to develop a new class of statistics 
to characterize the diffusion PDF. 

2.3. The orientation distribution function. An important tool in under- 
standing HARDI data is the orientational distribution function (ODF). The 
ODF quantifies the directional structure of the diffusion PDF in space. 
A popular object of study, it corresponds to several different functions in 
the literature. Tuch (2004) and Hess et al. (2006); Descoteaux et al. (2007) 
define the ODF to be 

i r°° 

(2.9) ODF T (0,0) = -/ a(ru)dr, 

Z Jo 

where x = ru, ||u|| = 1 and Z is a normalizing constant. Because this is not 
a true marginalization of a PDF (the increment needs a weighting by r 2 ), 
and weights lower scales heavily, the diffuse directional structure of the 
large-scale structure smooths the marginal PDF of orientations, giving it 
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Fig. 2. Diffusion processes displayed in both spatial and frequency (q-space) domains, 
with coloring representing density on the sphere. Ellipsoid diffusions are represented by 
their covariance matrix eigenvalues {A;}f =1 which govern a symmetric spatial decay, (a, g) 
Prolate (ellipsoid) diffusion process (Ai 3> A2 = A3). A prolate diffusion process is domi- 
nated by a single direction, represented by (a) a single peak in the diffusion PDF and (g) 
a great circle perpendicular to the diffusion direction in q-space. (b, h) Scalene (ellipsoid) 
diffusion process (Ai ~ A2 ^> A3). A scalene diffusion process has two competing directions, 
which makes the minor axes unequally matched in both spaces, (c, i) A mixture of prolate 
(ellipsoid) diffusions. This cannot be represented by a single unimodal diffusion PDF but 
must be represented by two directions, (d, j) and (e, k) These are both (nonellipsoid appear- 
ing) diffusion PDFs with asymmetric structure, suitable to model precursors to branching 
or forking (see text). Neither of these diffusion PDFs can be thought of as ellipsoid, (f, 1) 
Isotropic diffusion with no directional structure in space or q-space. 



a "blunted" appearance. A nonlinear transformation is necessary for the 
ODF to have a more peaked and clear directional structure. Wedeen et al. 
(2005) define the ODF as the truly marginalized PDF over all spatial radii 



(2.10) 



/•oo 

ODF W (0,0)= / r 2 a(ru)dr. 
Jo 
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An alternative version may be found in Jansons and Alexander (2003), where 
the orientational structure associated with a single radius is fitted to the ob- 
served data, that is, the persistent angular structure (PAS-MRI) algorithm. 
It is useful to note that the observed data are not associated purely with 
a single radius, and for this to be a mathematically correct procedure the 
observed HARDI measurements should be convolved with a suitable kernel 
prior to estimation. Despite this fact, the PAS-MRI method usually pro- 
duces good results in practice. All three of these orientational summaries 
are measuring different properties of the directional structure of the data, 
and only ODF\y(-, •) is a true marginal PDF. 

Another directional representation of diffusion data corresponds to the 
spherical convolution model [Tournier et al. (2004)]. In this model, q-space 
observations are modeled as convolved fiber ODFs, and fiber populations are 
estimated using deconvolution methods. The magnitudes are not comparable 
with previously-defined estimators of ODFs. Extensions to these methods 
have also been proposed: by modeling the ODF as a mixture of Bingham 
distributions [Kaden, Knosche and Anwander (2007)], and by regularizing 
the deconvolution problem by applying constrained optimization methods 
[Jian and Vemuri (2007)]. The solution in Kaden, Knosche and Anwander 
(2007) is parametric and the theoretical assumptions necessary to apply the 
regularized methods are, in general, violated [Jian and Vemuri (2007)]. 

The ellipsoid diffusion model (2.5) may be extended into a larger class of 
arbitrarily peaked and deformed diffusion PDFs by taking 

(2.11) A(x) = diag(Ai 1 (x),A 2 2(x),A 3 3(x)), A#(x) > OVx, 
with C a normalizing constant, to produce the diffusion PDF 

(2.12) a DE (x) = CVlAtYx^aiAtYxr^YxH), 

(2.13) a DE (T T x) = C^M^mK^r 1 / 2 ^). 

Because A(x) is a diagonal matrix, aD E (Y T x) exhibits the axes (1,0,0), 
(0,1,0) and (0,0,1). Applying a Fourier transform directly, with a change 
of variables, we note that the Fourier transform is mixed over the strengths 
in A(x), but exhibits the same orientational axes if the ordering in magni- 
tude of the eigenvalues does not switch over x. We have the model of 

(2.14) -4 DE (q) = C Iff lACxJp/afcdlACxJ-VaxIDe-^^^^x. 

This function can take the appearance of a deformed ellipsoid in space, 
and may then exhibit a different pattern of decay to the left and right of 
the dominant great circle in (/-space. For the regular ellipsoid distribution 
a E (x) if one eigenvalue is larger than the two others (say, Ai > A 2 > A3), 
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then the ellipsoid density [or equally in the case of the deformed density if 
inf x Ai(x) > sup x A2(x)] will observe a maximum at the values 



(2.15) q(/3) 



(3v 2 + yfT^v 3 , if /? e [-1,1], 

sgnQ3)(2 - |/3|)i, 2 

(2 - \P\) 2 v 3 , if/3G[-2,-l]U[l,2]. 



Figure 2g and h help to illustrate the behavior of (2.15), where the location 
on the "belt" is given by the value of f3. Note, the diffusion PDFs have 
been rotated in space compared to each other for a better visual perspec- 
tive. The maximum great circle in (/-space corresponds to the perpendicular 
vector ±i?i in space, where the diffusion PDF exhibits a maximum. The 
structure near the peak (x = ±v\) is mapped to a structure contiguous to 
the great circle, that is, q~ q(/3). Comparing the unimodal diffusion models 
(in Figure 2a, b and d), the microstructure of the diffusion PDF is mapped 
into behavior near or on the belt q(/3); see Figure 2g, h and j. The sca- 
lene structure of the diffusion PDF corresponds to variation on the belt 
(Figure 2h), while the asymmetry of Figure 2d and e are mapped onto the 
local structure of the delineation of the belt in Figure 2j and k. This mo- 
tivates us to investigate the structure of the diffusion PDF near the great 
circle of points {q(/3)} using distances from the great circle to characterize 
structure in the decay from the main peak. To obtain consistency in nota- 
tion, we define the set of points, or the great circle perpendicular to v, via 
G(v) = {q: f T q = 0, ||q|| = 1} and Q{v\) = {q(/3)}. It is convenient to keep 
both sets of notation for ease of exposition in the future. 



3. Scalar summaries and test statistics. 



3.1. Axes of symmetry. Before we can define appropriate scalar sum- 
maries in g-space, additional axes to the j3 axis (2.15) are required. For any 
fixed vector q(/3) £ Q{v{) we traverse a great circle using the vectors 

(3.1) qj_(a, (3) = avi ±y/l- a 2 q(/5), oG[-l,l], 

where for a G [— 2, 2]\[— 1, 1], the corresponding expression may be formed 
as in (2.15). Such a great circle for a fixed value of j3 will be referred to as 
a perpendicular great circle. 

An important component in the definition of our nonparametric sum- 
maries is the dominant great circle £(x max ) with x max given by 

(3.2) x max = argmax<^ <* A(q)dq\. 

v UqeGiy) J 

If «4(q) is an isotropic diffusion process, then x max is any vector in R 3 with 
a fixed norm. Alternatively, if «4.(q) is ellipsoid with Ai > A2 > A3, then 
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x max = v\. If there are two fibers, with relative weights of a\ and 0,2 of fiber 
populations with individual eigenvalues A*- 1 ) and A^ 2 \ then 




For example, if a± S> 02, then x max ~ v\ , or if a± = (12 = 1/2 and the great 
circles do not separate, then x max will lie precisely between the two max- 
ima of the two diffusion PDFs. Once the great circles start to separate the 
maximum will go with one of the two. 

3.2. Degree of nonuniformity. We represent a unidirectional Gaussian 
diffusion by plotting the value of A(q(fl)) (solid line) for [3 G [—2,2] in Fig- 
ure 3a. The magnitude on the dominant great circle is constant over dif- 
ferent values of f3 since A2 = A3. To illustrate the difference in variation 
across the dominant and perpendicular great circles, we also plot the value 
of A(q±(a, f3)) as a function of a for a fixed j3 (dotted line). This line per- 
fectly overlaps „4(q(/3)) at two locations, as it collides with the dominant 
great circle when it wraps around the sphere, and decays symmetrically 
from q(/3). 

We define a new coordinate system (a,/3), where we expect consistent 
variability in a and /3, using our parameterization of great circles (3.1). We 
plot the unidirectional Gaussian diffusion A(c[±(a, (3)) for all perpendicular 
great circles in the plane (Figure 3b). This prolate diffusion exhibits variation 
only in a, which is variation perpendicular to the dominant great circle. 
For the prolate diffusion example we can therefore reduce the variance by 
averaging across (5 and by considering the function strictly in terms of a. 
For comparison with the (a,/3) plane, the spherical representation of this 
Gaussian diffusion process is provided in Figure 3c. 

We use the one-dimensional great-circle summaries for a mixture of two 
Gaussian diffusions in Figure 3d, where the dominant great circle exhibits 
a large dynamic range relative to the perpendicular great circles. In fact, 
one can determine the number of peaks of the diffusion PDF by comparing 
the dynamic range of the diffusion between the dominant and perpendicu- 
lar great circles. For a complete picture we also represent the multi-modal 
diffusion in the (ct,f3) plane in Figure 3e, where variation is appreciable in 
both the a and (3 axes, and on the sphere (Figure 3f). 

To overcome the need to compare the variation along the dominant great 
circle with all perpendicular great circles individually, we define the average 
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(a) (b) (c) 




-2-1012 -2-101 



aorp P 

Fig. 3. One- and two-dimensional summaries of Gaussian diffusion processes in q-space, 
mapped onto the a and f3 axes (3.1) and their spherical representation, (a, b, c) Prolate 
diffusion process — eigenvalues (Ai S> A2 = A3), (d, e, f) Mixture of two prolate diffusion 
processes. The dominant great circle is the solid line in the one- dimensional summaries 
(a and d), while the dotted line is the diffusion from a single perpendicular great circle 
for (a) and the average perpendicular diffusion for (d). In the two-dimensional summaries 
(b and e) all great circles perpendicular to the dominant great circle are plotted on the y-axis 
to form the (a, /?) plane, and the final plots in (c) and (f) show the spherical representation 
on a single shell in Fourier space, corresponding to a fixed wave number magnitude. 



perpendicular diffusion via 

(3.4) A ± (a) = ^J^A(q ± (a,m))d$, 

with /3(t?) = cos(t9) for $ G [0,7r] and = — cos(-#) — 2sgn[cos(#)] denning 
/3($) for $ 6 [0, 2tt}. One may also define the average perpendicular diffusion 
over a half circle by prespecifying a fixed location on the dominant great 
circle and integrating in a window size ±1 around this location. This will 
prevent certain features being masked by the Hermitian symmetry of the 
(/-space measurements. If ^4(q) satisfies (2.2), then we have 

A±(a) = B((X ia 2 + (1 - a 2 )[X 2 \\q(P^)) T v 2 \\ 2 
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(3.5) 



+ A 3 ||q(/3(tf))^3|rD 



d-d. 



Thus, we are averaging the density function over small circles parallel to the 
dominant great circle and A±(a) measures the average diffusion at a given 
value of a. In the special case of A2 = A3, then 

1 r 2n 

(3.6) 



(3.7) 



A±(a) = — f B( v / \ 1 a 2 + \ 2 [l-a 2 ])d'& 
27r Jo 

= B(^X 1 a 2 + X 2 [l-a 2 }). 



The average perpendicular diffusion A±(a) provides a useful summary of 
variation perpendicular to the dominant great circle. We define a summary 
of the diffusion PDF via 



(3.8) 



max Q {^l_|_(a)} 



min a {^4j_(a)} 



max /? {.4(q 1 (0,/3))} 
min /3 {^(q ± (0,/3))} 



1. 



If the diffusion is isotropic, we know that Ai = A2 = A3. In this case we 
have -4_L(a max ) = A±(a min ) = B{y/\l) and >t(q(0, /3 max )) = -4.(q(0, /3 min )) = 
B(y/Xi), resulting in r = 0. If the diffusion is ellipsoidal and A2 = A3, then 
r = B{y/\2}/B{y/\]) — 1 > 0. If we adopt the mixture model, with multiple 
peaks, then it is possible to get t>0 even if we do not have a single 
diffusion PDF and we define 



(3.9) 



■4(q± («!,/?)) 

p ~ aT.aa i ^4(q± (a 2 , (3) ) 



mm max 



■A(q±(0,/W)) 



1. 



We note that under isotropy f = 0, while if we have a single ellipsoid diffusion 
f = r > 0. For a double tensor model f is more robust and will (in general) 
take on a lower value compared with r. In contrast to r and f, we could 
also study the variability in the g-space density directly in terms of the 
ODF. Tuch (2004), for example, defines the generalized fractional anisotropy 
(GFA) via 



(3.10) 



GFA 



™£r=i(ODF w (^)-l/n) 2 



in 



l)ELiODF 



w( 



i) 



1/2 



and this measures the nonuniformity of the spatial distribution, as do also 
the normalized entropy and the nematic order parameter [Tuch (2004)]. 
While the GFA quantifies the lack of uniformity in the ODF, if there is more 
than one fiber, determining its statistical properties is nontrivial, unlike the 
case for r and f. Another such measure, generalized anisotropy is defined 
in terms of the generalized trace of the tensor representation of the mean 
diffusivity [Ozarslan, Vemuri and Mareci (2005)]. 
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3.3. Measures of anisotropy. To determine the importance of the iden- 
tified dominant great circle (or orientation), we can, with a model of (2.2), 
compare B(y/Xi) to B^-y/X^) and B(y/X^). We define the following anisotropy 
statistic to perform such a comparison: 

log[AL(0)] _ Iog[S(v^)] 



(3 ' 11} ^ log[AL(l)] log[B(V*~i)Y 

where the last equality follows if A3 = A2. This statistic measures the degree 
of anisotropy over the g-space shell by comparing the peak-to-trough values 
(i.e., the value at the maximum great circle, compared to the value at the 
single point perpendicular to that maximum). Figure 3a displays the dif- 
ference between the maximum and minimum for an average perpendicular 
great circle. 

The decay ratio statistic quantifies the variability of the diffusion over the 
dominant great circle 

(3.12) C^max , '"^'f »' , . 

When the two smaller eigenvalues (A2 and A3) are approximately equal then 
£ ~ 1, otherwise £ 3> 1. The scalene diffusion in Figure 4c and d exhibits 
such structure (£ 3> 1). 

An indication of forking in white matter would correspond to an asymmet- 
ric decay of the diffusion PDF associated with different decays depending 
on the parity of the deviation. In this case we may no longer model the 
diffusion PDF as ellipsoidal. For example, in Figure 4a and b we see that 
while there is still a strong orientation from the dominant great circle, the 
PDF no longer exhibits symmetric decay away from the dominant great 
circle. Note that the decay is symmetric in a when averaged over the full 
sphere to produce A±(a). Hence, averaging over (3 is not appropriate if we 
want to detect asymmetry since a symmetric distribution will be obtained 
from the Hermitian symmetry of the HARDI measurements when averaging 
over a full great circle. A suitable asymmetry statistic to measure potential 
asymmetry is given by 

rm (l/2)/;/ 2 [^(q ± («W,/3))-.4(q ± (-«(^),/3))]^ 

(3.13) K{P) = T r. , 

(3.14) # max = argmaxK(/3(t?)), /3 max = /3(t? max ), 

(3.15) k= - / K{P(ti))dd. 

The definition of k is motivated by the wish to both obtain a test statistic 
with sufficient power and also to reduce its variance. The discrete approxi- 
mation to k will have a smaller variance than K{(3 max ). Asymmetry in the 
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(a) (b) (c) 




o or p P 



Fig. 4. One- and two-dimensional summaries of Gaussian diffusion processes in q-space, 
mapped out in the a and f3 axes (3.1) and their spherical representation, (a, b, c) An 
asymmetric diffusion process. This is apparent by the asymmetric decay in great circles 
perpendicular to the dominant great circle in the (a, (3) plane, (d, e, f) A scalene diffusion 
process with eigenvalues (Ai m A2 A3). The dominant great circle is the solid line in the 
one- dimensional summaries (a and d), while the dotted line is the the average perpendicular 
diffusion over j3 £ [—2,0] for (a) and all /3's for (d). The dashed line in (a) gives the 
average over all (3 's. In the two-dimensional summaries all great circles perpendicular to 
the dominant great circle are plotted on the y-axis to form the (ct,[3) plane. The final 
plots in (c) and (f) show the spherical representation on a single shell in Fourier space, 
corresponding to a fixed wave number magnitude. 

decay from the main peak may occur when the PDF is a mixture of diffusions 
with varying strengths. If the two populations are sufficiently separated and 
equivalent in magnitude, then this will be indicated by r and/or f and the 
diffusion will be recognized as a so-called "crossing fiber." If the mixture of 
diffusions contains two different strengths, then the dominating PDF will be 
recognized when determining x max . The remaining structure will not be fully 
consistent with a single tensor and will (in general) appear to be asymmetric 
compared to the dominant great circle. 

Let us discuss models that will lead to a different structure in the proposed 
summaries. We refer to Table 1 to summarize the properties of each statisti- 
cal test, and different diffusion PDFs lead to different structures. It may seem 
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Table 1 

The structure of the proposed diffusion models and the representation of their structure 
in terms of the proposed statistical summaries. Key to abbreviation where the statistics 
represent N-P/A (non-preference versus amsotropy) , C/E (circular versus ellipsoidal), 
S/A (symmetric versus asymmetric), I/M (isotropic versus multi-modal), M/U 
(multi-modal versus unimodal) 



Hypothesis 


Statistic 


Isotropic 


Prolate 


Scalene 


Mixture 


Heterogeneous 


N-P/A 


T 


small 


large 


large 


small 


large 


M/U 


f 


small 


large 


large 


small 


large 


I/M 




one 


small 


small 


small 


small 


C/E 


C 




one 


large 




large 


S/A 


H 




zero 


zero 




large 



insufficient to consider only an isotropic PDF, a single peak, a double peak, 
or something more heterogeneous. However, even with a fully parametric 
model of a Gaussian DTI framework, a two-tensor model has 13 (identifi- 
able) parameters and a three-tensor model has 19. If one considers acquiring 
60 gradient encoding directions (a common sample size), then one is forced 
to fit a highly-saturated model that results in noisy estimates — especially at 
higher 6-values where the orientational heterogeneity can be well resolved. 
Pushing much beyond a small number of parameters or features of interest 
is not advisable with such sampling. 

4. Estimation. 

4.1. Parameterizing the (a, j3) axes. Having proposed various summaries 
of the population of diffusion PDFs at a particular voxel, these must now 
be estimated from a set of diffusion measurements. The dominant direction 
may be estimated via 

t>l = x max = arg max <^ / A(q) dq \ 
u,IMI=i Uqea(u) J 

(4-1) 

= arg max FRT{^4} (x) , 

X 

where FRT{-} denotes the Funk-Radon Transform (FRT) as utilized in 
Tuch (2004). Note that «4.(q) refers to the multiresolution-based estimator 
[Olhede and Whitcher (2008a) and may be replaced by another appropriate 
estimator. We assume the availability of the quantity (<7*) 2 , an estimator of 
the variance of the error in -A(qfc) which we define to be a 2 . The variance of 
A(<ik) is assumed to be a 2 < a 2 and the variance of an interpolated value 
of the diffusion PDF is a 2 < a 2 < a 2 . The integral may be approximated 
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numerically by interpolating the observed HARDI measurements at evenly- 
spaced points along several great circles, each perpendicular to a given Xj. 

The effects of using different numerical methods for this step is a trade-off 
between increasing numerical accuracy and decreasing variance. Interpolat- 
ing using spherical harmonics reduces variance but may smooth out details 
depending on the choice of regularization; see the discussion in Descoteaux 
et al. (2007) and Hess et al. (2006). We instead use a locally linear inter- 
polation on the polar representation of the observed data that enforces the 
periodicity of the space in which it was sampled. Simple structures in terms 
of the observed points can mix over several spherical harmonics, and so the 
magnitude of individual spherical harmonic coefficients may not be large, 
even if the local coefficient is large. This fact makes the representation inap- 
propriate for using the smoothing methods of previous authors. The choice 
of interpolation procedure should be considered in terms of which statistic 
one is using, as the variance and bias must be balanced specifically for this 
purpose. We also note that spherical harmonics do not possess the same 
properties as Fourier vectors, and that an infinite number of harmonics are 
required for perfect reconstruction of a surface on a sphere, and so any 
reconstruction from the continuous basis will be inaccurate. 

The spatial maximum is determined from {FRT{„4}(xj)}j. The spatial 
location x max is an estimator of v\ and we estimate a vector in the linear 
subspace spanned by t>2 and V3 from the eigensystem of I — x max x^ ax , 
this yielding t>2 and 63, that maximize the difference in decay in the two 
axes. For numerical implementation we sample the estimated dominant great 
circle by discretizing a to {aj}j =1 and /3 to {Pk\k=l> ^ or an even integer N . 
A discretized version of (2.15) is then given by 



r (2 - p k )v2 



(4.2) q k ={ 



where 



y/1 - (2 - /3 fc ) 2 -u 3 , for k = -N/4, . . . , -1, 



PkV2 + yj 1 - $63, for k = 0,...,N/2-l, 

("2 - fa)V2 

- v / l-(2 + /3 fc ) 2 -u 3 , for k = N/2, 3N/4 - 1, 



2-cos(2tt£;/A0, for k = -N/4, . . . ,-1, 
(4.3) p k = { cos(2irk/N), for k = 0, . . . ,N/2 - 1, 

-2 - cos(2irk/N), for k = N/2, . . . ,3iV/4 - 1. 

We define q/% for the values of k not between k = —N/4, . . . , 3N/4 — 1 by 
cyclically extending (4.2). The choice of discretization guarantees the dis- 
tance between the orientation associated with a great circle and the great 
circle is one. The parameters aj and /3j are individually discretized to force 
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equal length increments on the great circle. Thus, we discretize the esti- 
mated dominant great circle via {qfc}. Once x max has been determined, the 
diffusion may be characterized directly in (/-space. Let us define the sampled 
great circle vectors for {qfc} in (4.2) via 



(4-4) qoj« 



where 



(2 - otj)vi 



v /l-(2-a j ) 2 qfc, 



atjij! + yj 1 - 
(-2-a 3 -)t>i 



for j 
for j 



-N/4, 
0,...,N/2 



Vl-(2 + a,) 2 qfc, for j = N/2, 3N/4 - 1; 



(4.5) 



a : 



2-cos(2nj/N), 
cos(2vr j/N), 
-2-cos(2irj/N), 



for j = -N/4, . . . , -1, 
for j = 0,...,N/2-l, 
for j = N/2, 3N/4 



1. 



A discretized version of the average perpendicular diffusion (3.4) is given 
by A±(aj) = N~ 1 ^2 k A(<i_Ljk), and we can sum over any N consecutive k 
(e.g., it does not matter exactly how we sum over k) because of the periodic 
extension. 



4.2. Diagnosing nonuniformity. In order to test large-scale properties of 
the diffusion directly in (/-space, we consider the following statistical hypoth- 
esis Hq :^4(q) = A Vq versus Hi :A(q) = «Ae(<1)- Our test statistic is based 
on a discretized version of (3.8), given by 



(4.6) 



max.j{A±(otj)} 
minj{„4._|_(aj)} 



maxfcj^l(qfc)} 
mm k {A(q k )} 



1. 



The distribution of this test statistic is derived in Supplementary Material. If 
the observations are isotropic, then the properties along the dominant great 
circle will be equivalent to the properties on the perpendicular great circle 
(ignoring any random/discretization errors). The estimators of A and a, 
under the null of -4(qfc) — A + as, are given by 



(4.7) 



(4.8) 



N 1% 27r Jg&i) 

a A = VpMAD{l(qfc) -A N :k = l,...,N}, 



where < p < 3, qfc defined by (4.2) and MAD{-} is the maximum ab- 
solute deviation. These equations provide estimators of the mean value of 
the isotropic diffusion and the standard deviation of -4(q) at the observed 
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measurements. The parameter p is a constant depending on the linear in- 
terpolation method used for the implementation. Taking a value of p = 3 is 
suitable for our choice of numerical interpolation and we define 

(4.9) U = t4^. 

We can compute the critical value u a using the fact that Fu(u a ) = 1 — a, 
where Fjj(-) is defined in Supplementary Material. We report two critical 
values here, «o.05 = 0.1185 for the m which is consistent with our sampling 
scheme, and UqqJ^ = 1.9637 with a conservative distribution approximation. 

We also develop a new test based on a null of a multi-modal diffusion, 
where we define multi-modal in terms of (^4 m ax-^min)/(-4min-4max) < c = 2, 
where A mSiX and A m \ n are the maximum and minimum on the perpendicular 
great circle minimizing (3.9) in f3, respectively, while A ma , x and -4 m i n are the 
maximum and minimum on the dominant great circle. The value of c is 
fairly arbitrary, but to develop a powerful method of separating the clearly 
unimodal from the multi-modal, some level must be chosen based on the 
deterministic structure of the sampled diffusion PDF. To separate unimodal 
from multi-modal PDFs, we start from f and define 



(4.10) T = min< max 



±hk) 



lnd2 M(q ±i2fc ) 



max fc {^l(qfc)} 
min fc {^l(q fc )} 



1. 



We shall now choose to distinguish the multi-modal from the unimodal, 
and so normalize the statistic using U = (T — [c — l])A m i n / (<r_4\/2c 2 + 2), 
where A m i n = A(vi). The distribution of this test statistic is derived in 
Supplementary Material, under the specified null hypothesis. 

If, on the other hand, we have failed to reject the null hypothesis "A(q±(a, 
(3)) equally variable in /3 for a = as it is in a for a fixed /?," then based 
on the T-statistic we need to distinguish voxels that indicate two fiber pop- 
ulations versus isotropic voxels. Let us define a discrete version of (3.11) to 
be 

(4.11) x= l ^4M. 

iog[.4 ± (i)] 

We can interpret £, and the sample version X , as the degree of anisotropy 
from the average perpendicular great circle. We recognize that the statis- 
tic X is comparing the average apparent diffusion coefficient (ADC) on 
the great circle to the average ADC perpendicular to the great circle, or 
that (4.11) may be rewritten in terms of the ADC at a fixed value of b via 

(4.12) X = C(q k ) I £(4uv/4fc)- 
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The ADC is C(cjj) = — b^ 1 log A(c\j), for a defined set of q vectors q,- [Alexan- 
der, Barker and Arridge (2002)]. With an assumption of ellipsoidal struc- 
ture (2.2) we have averaged the ADC to reduce variance when estimating £ 
without accruing bias. We define Xk as the sample anisotropy calculated 
using only the kth. perpendicular great circle via 

(4.13) X k - l0g[I((lfc)] 



log[^(q±iv/4ifc)] 

and refer to (4.4) for the definition of q_i_jfc- Under moderate-to- high SNR 
may be approximated by a Gaussian random variable. We quantify uncer- 
tainty, when there are potentially several peaks, using <5"2 = min{<74, <r*}, 
where oj^ is defined in (4.7) and a* is the available estimator for a. By 
using the minimum, we ensure that the estimated variance is not inflated 
compared to its pre-smoothing value. 

For those voxels where isotropy cannot be rejected, we may now distin- 
guish between isotropy and a multiple-tensor model using X. The multi- 
modality statistic is given by 

(4.14) Q= p ( X ~ 1 " > \A N lagA N \. 

So we consider the test Ho\A(q) =A Vq versus H\ : maxfc{.A(qfc)} 3> 
minfc{^l(qfc)} (multiple peaks) and use Q as the test statistic, whose dis- 
tribution under the null is provided in Supplementary Material. The three 
tests outlined here allow one to at a single voxel diagnose the structure 
of the diffusion PDF, where U is used to separate anisotropic PDFs from 
isotropic PDFs, U is used to separate ellipsoid PDFs from multi-modal PDFs 
and Q is used to separate multi-modal PDFs from isotropic PDFs. 

4.3. Diagnosing asymmetry. Having established methodology to discrim- 
inate the number of peaks in the diffusion PDF at a single voxel, we now 
provide additional methodology to characterize the diffusion PDF as scalene 
versus other forms of asymmetry, for example, to observe the indication of 
forking or fanning white-matter structure. Let us define 

/„ 1k n , log-4(qfc) 

(4.15) fcmax = arg max 



l<k<N/A log^(q fc+w/4 ) 

The effective degrees of freedom parameters (m,m') are related via m < 
w! < 2m, for robustness, so that 

^ 416 ^ z _ lo g^(qfc m ax+7V/(2m')) 



log^(q fcmax+Ar/ ( 2 . m /) +7V/4 ) 
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We remark that Z is related in some sense to Xk (refer to Figures 3 and 4) . 
The statistic Xk compares the value of the diffusion PDF at location k 
on the dominant great circle (a = 0) to the value at the perpendicular to 
the dominant great circle («/0). The statistic Z in contrast looks at the 
difference in values of the diffusion PDF on the great circle itself (a = 
and /3 varies). Under the null hypothesis A(-) is constant on the great circle, 
and if the medium and minor eigenvalues are approximately equal, then 
E{Z} = £ ~ 1, otherwise £ 3> 1. A normalized version of the decay ratio 
statistic (4.16) is given by 

(4.17) v= (Z-l)\A N logA N \ ^ 

<5"2 

A suitable threshold for this statistic may be found in Supplementary Material. 
The statistics Q and V, used to test different hypotheses of nonisotropic de- 
cay, have similar forms. 

The summary statistic k in (3.15) allows one to diagnose white-matter 
microstructure that is not consistent with a single ellipsoid diffusion. De- 
partures from such a single ellipsoid diffusion model may be attributed to 
partial- volume effects, or a heterogeneous population of fibers [Behrens et al. 
(2007)]. For such a model (2.2) is no longer appropriate and we would rather 
fit a mixture model with unequal populations — or possibly -4de(0- In such 
circumstances one cannot use the average perpendicular great circle to un- 
cover asymmetry since averaging over all possible /3's will produce a sym- 
metric distribution regardless of the underlying fiber characteristics. Taking 
fc max = argmaxfc{Pfc}, we define the asymmetry statistic via 

fcmax + A78 

< 418 » K = mr\ E ^ 

fc = fcmax-A r /8 



£^=i-4(q±jfc) 



(4.19) P k 



Full details on the distribution of K may be found in Supplementary Material. 
We have chosen N/8 to improve the power, since averaging decreases the 
variance, but the asymmetry is greatest near the maximum (compare with 
Figure 4b). For tests at a specific voxel we perform the hypothesis test 
Hq : k = versus H\ : k ^ 0, using quantiles from the standard Gaussian PDF 
</>(•). This text identifies diffusion PDFs that are non-Gaussian in terms of 
the parity structure in the principal axes. However, it does not compare 
the maximum and minimum of a perpendicular great circle, rather it finds 
a set of perpendicular great circles for which the decay around the domi- 
nant great circle is asymmetric and estimates this average asymmetry, for 
example, Figure 4b. 
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The usage of the statistics is now combined at a voxel level. The most 
important step is to classify the voxel as isotropic, unimodal or multi-modal. 
With this information the local structure of the peaks in the diffusion PDF 
may be further characterized. This is similar to the situation when com- 
paring PDFs between voxels for fiber tracking (tractography) where the 
number of mixture components in the voxels is the first priority and then 
the components in the diffusion PDF are matched between voxels using local 
characteristics corresponding to structures at fine scales. From knowledge of 
white-matter structure one would anticipate varying values of asymmetry 
before a forking fiber structure, and this would allow us to smoothly go be- 
tween single voxels with unimodal diffusion to mixtures. These topics shall 
be discussed in subsequent sections. 

4.4. Example: Crossing and forking fibers. We consider two typical het- 
erogeneous white-matter structures, a crossing fiber and a forking fiber (Fig- 
ure 5). The spatial representation of the forking fiber is provided in the first 
row of the figure, denoted by (i), while the g-space representation of the 
forking fiber is given on the second row of the figure, denoted by (ii). In the 
spatial representation we see a single fiber population in voxel (i,a) and, as 

(a) (b) (c) (d) (e) (f) (g) 



%) %) d 9 



(") 



(i") 



<• * 3 w * * 



(iv) 



Fig. 5. An illustration of the evolution of a diffusion PDF through a number of adjacent 
voxels in space and q-space. The first and second rows are the spatial and q-space evolution, 
respectively, of the diffusion PDF for a forking fiber. The third and fourth rows are the 
spatial and q-space evolution, respectively, of the diffusion PDF for a crossing fiber. The 
aim of the plots is to show the changing q-space structure over this evolution. 
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we traverse from left-to-right, the two populations become more apparent 
until a second fiber appears in voxel (i,g). The g-space version of these two 
populations shows a scalene distribution developing in voxels (ii,b)-(ii,e). As 
the forking progresses, from left-to-right, it appears both highly warped and 
scalene until the distribution clearly displays multiple fibers in voxel (iii,g). 
These fibers were generated by aggregating two densities via 

i(q 1 () = fl 1 (tM 9 (q;A( 1 )(f),TW(f)) 

(4.20) 

+ (l-a 1 (t)M 9 (q;A( 2 )(t),T( 2 )(t)). 

In the case of the forking fiber, a±(t) = 1 — i/2 and the main directions of 
Y«(t) and r&(t) are given by (1,0,0) and (cos(7rt/2), sin(7rt/2), 0), re- 
spectively. The individual tensors take values similar to Ai(q) and t £ [0, 1]. 
The spatial representation of the crossing fiber is the third row (iii) of Fig- 
ure 5, with its corresponding q-space representation in the fourth row (iv). 
The ellipsoid appears prolate in voxel (iii, a), then two fiber populations 
are present in voxel (iii,d) and eventually the fiber population returns to 
a prolate shape. With respect to the parameterization of the crossing fiber 
in (4.20), ai(t) E {1,0.75,0.5} and the two fibers cross at 90 degrees with 
parameters similar to „4i(q). 

The description of a crossing fiber is in many ways simpler than a forking 
fiber. Table 2 lists the summary statistics (r,f, k) for the crossing and 

Table 2 

Discretized summaries based on nonparametnc measures of symmetry for a modeled 
forking and crossing fiber, compare with Table 1. The summaries show how the statistics 
evolve over a sequence of voxels undergoing forking or crossing 











Forking fiber 








Statistic 
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(i,c) 
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(i,e) 


(i,f) 


(i,g) 


T 


10.18 


5.12 


3.35 


1.86 
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0.07 
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8.98 


4.94 


3.36 


2.13 


1.19 


0.33 


-0.07 




0.12 


0.18 
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0.27 


0.35 


0.53 


0.69 
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1.03 


1.51 


1.72 


1.91 


2.06 


2.21 


2.67 
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-0.01 


0.03 


0.17 


0.35 


0.45 


0.54 


0.30 
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9.19 
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1.26 


-0.07 
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0.69 


0.32 


0.12 
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1.04 


1.04 


1.74 


2.67 


1.74 


1.04 
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0.00 


0.00 


0.09 


0.30 


0.09 


0.00 
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forking fiber examples in Figure 5. Note that these deterministic summaries 
have not been normalized, unlike the statistics in Section 4 (since there is 
no noise with which to compare). The mixture of populations of unequal 
strength in the forking fiber shows a number of characteristics not found 
in the crossing fiber. For example, the forking fiber is clearly diagnosed as 
a single population until voxels (i,e)-(i,g), where there is increasing hetero- 
geneity in the fiber population. This is exhibited by increasing values for 
the decay ratio ^-statistic, and the asymmetry K-statistic. For the crossing 
fiber, we clearly detect the multiple-fiber population in voxel (ih,d) using 
either the r or f statistics. The multiple-fiber characteristics in the forking 
example are more complex, where the second fiber population is initially 
dominated by the first. If we examine the crossing fiber more closely, there 
is little apparent asymmetry and we can compare the asymmetry statistic, 
where k ~ versus 0.15 < k < 0.45 for voxels (i,c)-(i,e). To distinguish mul- 
tiple fibers from uniformity, we observe that £ < 1, which is the expected 
value under the hypothesis of isotropy. 

5. Simulation study. We illustrate the properties of the proposed g-space 
summary statistics for the diffusion PDF on a variety of simulated diffusions 
processes. The following six models attempt to cover common, and not so 
common, diffusion processes that include both unimodal and multiple ten- 
sors: 

A(q) =exp(-tq T Diq), i = 1,2,3, 
Di = 68eief + 8e 2 &f + 8e 3 eT, 

(5.1) 

D 2 = 68§ief + 15e 2 e 2 + e 3 ej , 
D 3 = 28eief + 28e 2 ef + 28e s el, 

Ai(q) = 0.5exp(-tq T Diq) + 0.5exp(-iq T D 4 q), 

(5.2) 

D 4 = 68e 2 e 2 + 8eief + 8e 3 e 3 , 

A(q)=exp(-m|q T e 2 | 2 ) 
(5.3) x |exp(-68t|q T ei| 2 ) x [exp(-0.2t|q T e 3 | 2 ) + exp(-35t|q T e 3 | 2 )] 

+ -,D(v / 68tq T e 1 )[D(v / 35tq T e 3 ) - D(V0Mq T e 3 )}\, 

7T 

A(q) = 0.3exp(-tq T Diq) + 0.7exp(-tq T D 5 q), 

(5.4) 

D 5 = 42.5ei&f + 14e 2 e 2 + 20e 3 e 3 , 
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where D(x) = exp(— x 2 ) Jq exp(i 2 ) dt is the Dawson function [Abramowitz 
and Stegun (1972)]. We define ej = IZej, where the matrix 1Z rotates the 
axes (ei,e2,e3) to a new coordinate system (ei?^,^). This extra step is 
added to protect against systematic bias in our estimation procedure due to 
the diffusion PDF coinciding with the sampling grid. In Aq(c[) this rotation is 
not implemented, but (61,62,63) has been rotated with respect to (ei, 62,63) 
to produce an asymmetric diffusion in the multi-tensor model. 

These diffusion processes are displayed in both spatial and frequency do- 
mains in Figure 2, where A±(-) is a prolate diffusion model (a and g), A2(-) 
is a scalene diffusion (b and h), An{-) is a mixture of two crossing fibers 
(c and i), Aq(-) is the first asymmetric diffusion (d and j), and A§(-) is 
the second asymmetric diffusion (e and k). For completeness an isotropic 
diffusion model As(-) is shown in Figure 2f and 1. 

We have chosen to define Dj = 4 x 10 10 Dj, i = 1, . . . ,4, and normalized 
||q|| =1. With t = 0.04 this corresponds to b = At x 10 10 = 1600 s/mm 2 and 
the trace of the first three nonnormalized matrices Dj as 2.1 x 10~ 9 m 2 /s 
[Alexander (2005)]. The function -4s(q) is obtained from the magnitude 
of the FT of an asymmetrically decaying diffusion process in space. We 
illustrate a range of behavior for the scalar statistics defined in g-space using 
these test functions, providing only a subset in order to compare and contrast 
their performance. We simulate 1000 realizations for each test function and 
add Gaussian noise with standard deviation of ^4(0) /2, A(0)/10, „4(0)/20 
and ^4(0)/30 to both the real and imaginary channels using a 60-direction 
HARDI sampling scheme. 

Results, provided in Table 3, are consistent with varying degrees of the 
SNR. The prolate diffusion A± is clearly detectable, down to an SNR = 1/10, 
despite using nonparametric methods via the [/-statistic. Detecting the sca- 
lene diffusion depends on the SNR, while the isotropic diffusion is clearly 
distinguishable from its alternatives under the full range of SNR using the 
[/-statistics. The multi-tensor diffusion A4 is difficult to classify using the 
[/-statistics and its correct classification depends on how well the location 
of the dominant peak is estimated. If the dominant peak is well determined, 
then the [/-statistic clearly recognizes the density as anisotropic, if not, the 
g-space measurements are characterized as non-Gaussian instead of multi- 
modal. If one was only concerned with empirically separating prolate diffu- 
sion PDFs from multi-modal diffusion PDFs, rather than performing a hy- 
pothesis test, then this would be relatively straightforward, for example, 
retaining 95% of the unimodal Gaussian with the SNR= 1/20 leads to re- 
jecting all but 11% of the multi-tensor realizations (see the [/-statistic). 
Since we are interested in detecting ellipsoidal decay around a single direc- 
tion, the variation over the dominant great circle will be large for anisotropic 
voxels with ellipsoidal decay as well as for multi- modal diffusion PDFs. At an 
SNR = 1 /20 the [/-statistic provides complimentary information by strongly 
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Table 3 

Hypothesis tests for the six diffusion processes {Ai}i=i, specified in (5.1), where the 
number of rejected hypothesis are provided as a single number out of 1000 tests, or as 
a fraction if fewer than 1000 tests were performed. The nominal size of the tests is 5% 
for U , V and Q, while the nominal size is 10% for K and U . The hypothesis tests have 
been carried out at different SNRs, where the noise standard deviation is increasing as 
you go further down the entries SNR 6 {1/30, 1/20, 1/10, 1/2}. Keys to the abbreviations 
are N-P/A (nonpreference versus anisotropy), C/E (circular versus ellipsoidal), S/A 
(symmetric versus asymmetric), I/M (isotropic versus multi-modal) and M/U 
(multi-modal versus unimodal) 



H /H 
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i /on 
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1000 


988 


26 


492 


1000 


802 


C/E 
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146/1000 


924/988 


0/26 


382/492 


120/1000 


495/802 


S/A 




K 


191/1000 


108/988 


10/26 


258/492 


491/1000 


208/802 


I/M 




Q 


0/0 


12/12 


21/974 


420/508 


0/0 


195/198 


M/U 
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991/1000 


136/988 


0/26 


239/492 


996/1000 


19/802 


SNR = 


1/20 
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1000 


855 


26 


484 


1000 


727 


C/E 
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153/1000 


794/855 


0/26 


338/484 


148/1000 


267/727 


S/A 
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148/1000 


38/855 


10/26 


199/484 


331/1000 


136/727 


I/M 
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0/0 


0/145 


23/974 


259/516 
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201/273 


M/U 







945/1000 


46/805 


0/23 


151/484 


942/1000 


8/727 
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1000 


239 
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441 


998 


457 
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225/1000 


214/239 


0/34 


192/441 


194/998 


58/457 


S/A 
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89/1000 


1/239 


10/34 


99/441 


139/998 


74/457 


I/M 
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449/761 


20/966 


20/539 


1/2 


18/543 
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239/1000 


7/239 


0/34 


27/441 
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4/457 
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45 


21 


25 


31 


37 


24 
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2/45 


2/21 


5/25 


6/31 


4/37 


4/24 
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5/45 


4/21 


2/25 


1/31 


4/37 


3/24 


I/M 
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1/955 


3/979 


2/975 


2/969 


4/963 


5/976 
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2/45 


0/21 


1/25 


0/31 


5/37 


1/24 



separating the prolate (94% rejected) from the multi-tensor model (15.1% 
rejected, near the nominal value of 10%), but fails to distinguish between the 
scalene and the multi-tensor models (Table 3). The highly scalene diffusion 
is mistaken (not surprisingly) for a multi-modal diffusion and such structure 
may be approximated using two tensors, especially when sparsely sampled 
on the sphere. 

The two distributions with constant behavior on the dominant great circle 
are not diagnosed with asymmetric decay, while the null hypothesis is re- 
jected for A2 in a substantial number of cases in Table 3. The misdiagnosed 
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multi-modal diffusion PDF A& also has the null hypothesis of multimodal- 
ity rejected for a substantial number of cases. This is to be expected since 
the observed diffusion will experience considerable variation over the domi- 
nant great circle, consistent with observing a diffusion process with a single 
dominant direction and ellipsoidal decay. 

We fail to reject the null hypothesis of symmetry for the two diffusion 
processes that are symmetric (A± and A2) in most cases, while we reject 
a larger proportion for ^.5. There is unfortunately a lack of power in this 
test which is due to sampling 60 directions, limiting the performance of the 
test statistic. For Aq and keeping the SNR = 1/20, we reject the null hy- 
pothesis 38.9% with 60 directions. For A5 we reject the null 51.6% of the 
time using 245 directions at SNR= 1/20 — a clear increase from 35.2% with 
60 directions. Increasing the SNR also increases our power to detect such 
asymmetry, as shown in Table 3. The power of the test improves as the num- 
ber of directions increase or the amount of asymmetry (better characterized 
with better spherical sampling) increases. For all of the structural tests per- 
formed here there is a direct similarity in effect of increasing the number of 
grid points to improve the size of the mean under the alternative hypothesis 
or directly decreasing the variance. This is because as the mean of the al- 
ternative hypothesis increases with improved sampling of g-space, this has 
the same effect as increasing the SNR, as the test statistics are (approxi- 
mately) functions of their ratio. This direct exchangeability of sampling in 
frequency versus SNR holds until the distributional approximations break 
down because of poor resolution in g-space or a diminished signal-to-noise 
ratio. 

6. Analysis of clinical data. HARDI data were acquired from one normal 
subject (30 year old, male Caucasian) in a Siemens TIM Trio 3.0 Tesla 
scanner using a 32-channel head coil. Measurement of 64 gradient directions 
(6 = 1600 s/mm 2 ) and one T2 image (b = 0) were obtained using a twice- 
refocused diffusion preparation. The slice prescription was 64 slices acquired 
in the AC-PC plane, TE = 95 ms, FoV = 240 x 240 mm, base resolution = 
128 x 128, slice thickness of 1.9 mm and cardiac gating was applied. 

Regions of interest (ROIs) from two slices of the clinical data are pro- 
vided to illustrate the statistical summaries developed in this paper. Slice 1 
contains an ROI that is dominated by single-fiber voxels containing struc- 
tures such as the corpus callosum and cingulum. Figure 6a shows the voxels 
using the common color-coding convention [i.e., RGB for the (x,y,z) coor- 
dinates] weighted by the estimated fractional anisotropy (FA) at each voxel. 
The FA for the ROI is reproduced in Figure 6b along with the p-values for 
the anisotropy and ellipsoidality statistics in Figure 6c and d, respectively. 
We select a very liberal threshold (p = 0.15) for the purpose of exploratory 
data analysis, not confirmatory data analysis. We observe very few voxels 
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Fig. 6. Axial slice from clinical HARDI acquisition. Color-coded fractional anisotropy 
(FA) for the whole slice is displayed in (a) along with the boundaries for the ROI (Re- 
gions of interest). For the zoomed-in ROI: color-coded FA (b), anisotropy p-values (c), 
ellipsoidality p-values (d), unimodality test statistic (e) and multimodality p-values (f). 
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that indicate asymmetry at specific voxels, while the ellipsoidality p-values 
indicate quite a few voxels that exhibit prolate diffusion. These voxels are 
located at the borders of strongly directional structures such as the corpus 
callosum and cingulum, and reaffirm the results obtained in the simulation 
studies. Additional information about the structure is obtained by plotting 
the test statistic for unimodality and the p- values from the multi-modality 
test statistic in Figure 6e and f, respectively. The corpus callosum, and to 
a lesser extent the cingulum, produce large values in the unimodality test 
statistic as to be expected from those structures. Multimodality is detected 
in voxels with reduced FA and/or on the edges of prominent white-matter 
structures. The pattern of multi- modal voxels identified in Figure 6f in gen- 
eral do not appear to overlap with those voxels that were identified using the 
ellipsoidality statistic, providing evidence that this methodology is detecting 
distinct features in the white-matter microstructure. 

The ROI selected in slice 2 captures more complicated interactions be- 
tween white-matter structures such as the corticopontine tract, anterior 
thalamic radiation and corpus callosum (Figure 7a). The FA for the ROI 
is reproduced in Figure 7b along with the p-values for the anisotropy and 
ellipsoidality statistics in Figure 7c and d, respectively. Asymmetry is dif- 
ficult to detect in these data, but ellipsoidality is quite apparent along the 
boundaries of the corpus callosum and around the projections into gray 
matter. The test statistic for unimodality in Figure 7e complement the el- 
lipsoidality results quite well, picking out dominant prolate diffusion (e.g., 
the voxels dominated by the corpus callosum and to a lesser degree the 
cingulum) around which the ellipsoidality measure is finding more complex 
voxels. Finally, the test statistic for multimodality in Figure 7f clearly iden- 
tifies voxels where the three dominant white-matter structures in this ROI 
converge, and all other statistics fail to detect any specific structure. The 
statistical summaries developed here provide complementary information 
about white-matter microstructure in clinically acquired data. 

We focus on a few specific voxels in Figure 7 using the Funk-Radon 
Transform (FRT) without smoothing. As recommended by Tuch (2004), 
we have taken the standardized raw FRT to the power five to emphasize 
structure in the display. Figure 8a and b show the two most anterior voxels 
that are plotted in Figure 7 (indicated by yellow dots) . This tract appears to 
be "bending" as we move from anterior to posterior, indicated by the shift in 
direction of the dominant direction seen in the FRTs. The statistics quantify 
this behavior; the p- values for asymmetry are 0.14 and 0.02 respectively 
(indicating that the posterior-most voxel is bending more). The unimodality 
of the anterior voxel is seen from the large unimodality statistic in Figure 7e. 
We then look at a voxel in a more heterogeneous area, where the major fiber 
tracts appear to merge: the statistics here indicate multi-modality dominates 
as is seen in Figure 8c and backed up by Figure 7c-f. We observe the most 
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Fig. 7. Axial slice from clinical HARDI acquisition. Color-coded fractional anisotropy 
(FA) for the whole slice is displayed in (a) along with the boundaries for the ROI. For the 
zoomed-in ROI: color-coded FA (b), anisotropy p-values (c), ellipsoidality p-values (d), 
unimodality test statistic (e) and multimodality p-values (f). 
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(a) (b) (c) (d) 

© G © O 

Fig. 8. The raw Funk Radon Transform (FRT) from a collection of voxels indicated 
by yellow dots in Figure 7. These are plotted in order of decreasing xi- coordinate (or 
going from the top of the image to the bottom). Subplots (a) and (b) both reject the null 
hypothesis of no asymmetry, with (a) not rejecting prolate diffusion in favor of scalene 
diffusion. Subplot (c) rejects unimodality in favor of multimodality and also rejects isotropy 
in favor of multimodality. Subplot (d) is unimodal. The raw FRTs are consistent with these 
diagnoses. 

central voxel has summary statistics that are ellipsoidal but not asymmetric, 
clearly observed in Figure 8d. The clinical data have provided evidence at 
a voxel level, backed up by statistical hypothesis testing and observed in 
the FRT visualizations, that interesting white- matter microstructure may 
be detected and characterized using the methodology proposed here. 

7. Discussion. We have introduced a new set of tools for characterizing 
orientational structure from HARDI measurements directly in g-space. This 
methodology is unique when compared with existing methods that rely on 
reconstructing the spatial information from (/-space by different methods 
of marginalizing the spatial distribution, that is, from calculating a spatial 
ODF. An ODF has a different meaning if calculated directly from a Gaussian 
model, from the nonparametric FRT (average orientational distribution over 
all radii without using the correct volume increment for a marginal PDF) or 
using PAS-MRI (orientational distribution associated with a single spatial 
radius or scale). In general, the magnitude associated with an ODF is not 
comparable between methods, neither is the distribution of noise artifacts. 
Our methodology is technically linked to the FRT, but unlike the FRT we are 
not constrained to scalar measures calculated from averages on great circles 
in g-space, and our methods do not depend on appropriate marginalization 
to produce summaries. The interpretation of our statistical summaries is 
straightforward, but we note that in improvements in data acquisition, such 
as increased sampling of directions. 

Most established methods for characterizing features in white-matter mi- 
crostructure have focused on the problem of determining the number and 
orientation of peaks in the diffusion PDF. None of the "magnitude" infor- 
mation of these solutions is comparable or indeed interpretable apart from 
DTI-based models. Savadjiev et al. (2006) have already commented on the 
unsuitability of such magnitudes as quantitative measures. The problem with 
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this fact, and the nonlinear transformation often employed for representing 
(/-ball estimates, is that the coherent treatment of noise artifacts becomes 
much more difficult. The advantage of our theoretical framework, developed 
for summary statistics, is that we may perform hypothesis tests using crit- 
ical values that are not functions of unknown parameters. We stress that 
simulation studies for features of diffusion PDFs are in general misleading 
unless the proposed summaries are true statistics, that is, their distributions 
under null hypotheses are parameter independent. For example, critical val- 
ues determined from Monte Carlo studies for a given diffusion PDF will 
not (in general) be applicable to other diffusion processes than the simu- 
lated process since these critical values are parameter dependent. This can 
be compared to calculating a simple mean rather than a t-statistic. If we 
try to elicit the distribution of the sample mean using simulations at fixed 
variances, then these critical values are only useful for variables with the 
same variance. 

Various nonparametric procedures have been proposed to summarize 
HARDI data using more than its estimated orientation, for example, by 
investigating the model order of spherical harmonic decomposition [Frank 
(2001); Alexander, Barker and Arridge (2002); Descoteaux et al. (2006)]. 
Chen et al. (2005) modeled the ADC using a product of a truncated spher- 
ical harmonics series. In general, an infinite order of spherical harmonic 
terms must be taken to approximate an arbitrary Gaussian mixture, but 
they argued that a crossing fiber should be sufficiently reproduced by such 
a truncated representation and expressed its complexity using the normal- 
ized terms in the spherical expansion. Other representations include ex- 
pressing the ADC in terms of higher-order tensors and spherical harmonics 
[Descoteaux et al. (2006)], or just via a spherical harmonic representation 
[Frank (2001)]. Second-order terms in a spherical harmonic decomposition 
contribute to describing a single-tensor fiber, but more complicated struc- 
ture must be described in terms of corresponding spatial properties of the 
PDF directly, rather than the fourth- and higher-order terms which give too 
much freedom in structure to be a precise tool for the description of fine 
spatial features. Other measures of the entropy of the diffusion PDF have 
been proposed by Rao et al. (2004). 

Rather than solely focusing only on the number of peaks in the diffu- 
sion PDF, we have characterized white- matter microstructure through the 
diffusion PDF directly in (/-space without parametric assumptions or impos- 
ing smoothness constraints, as we use a variable bandwidth estimator rather 
than employing a fixed bandwidth smoother [Olhede and Whitcher (2008a)] . 
The tissue microstructure is identified as variation in summary statistics 
that deviate from a simple, symmetric model for the diffusion PDF and is 
characterized in behavior relative to the identified dominant great circle in 
(/-space. Ellipsoidal diffusion PDFs (2.2) are simple in structure and imply 
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the existence of a dominant great circle. The deformed ellipsoid class is less 
stringent in structure, and permits asymmetric decay in minor axes — for 
example, (2.14) — while still conforming to the existence of a dominant great 
circle. We describe the precursor to forking structures by either a deformed 
ellipse or a mixture model, to capture further asymmetric structure. We dif- 
ferentiate between different white-matter microstructure by examining vari- 
ation over that great circle, or variation perpendicular to the great circle. 
Allowing for a greater variety of structure in a unidirectional diffusion PDF 
implies that the power to detect multi-modal diffusion is necessarily reduced 
compared to using a parametric multi-model model, if the proposed para- 
metric model is correct. We characterized single peak densities by additional 
summaries, such as the anisotropy statistic, the decay ratio statistic and the 
asymmetry statistic. The synthetic forking fiber in Figure 5 shows an evolu- 
tion of such measures as we go between a single fiber and a forking fiber. The 
synthetic crossing fiber in Figure 5 does not exhibit the same asymmetries. 

If one enforces a strict Gaussian (single diffusion tensor) model, then all 
variation away from symmetry around the dominant direction will be inter- 
preted as evidence for a multi- modal diffusion [Parker and Alexander (2005); 
Hosey, Williams and Ansorge (2005); Behrens et al. (2007)]. Modeling us- 
ing non-Gaussian PDFs allows us to fit asymmetric structure, rather than 
just the model indicating a lack of fit of a single peak. However, using such 
models leads to a loss of power if a Gaussian mixture model is appropriate. 
Caution should be exercised in order to protect against over-interpreting fit- 
ted models. With a model that only includes a family of mixtures of Gaus- 
sian diffusion processes, one is constrained to estimate a Gaussian mixture, 
however, for a small number of sampled directions there will inevitably be 
issues with identifiability. The same realizations may in some cases equiva- 
lently be derived from a unimodal diffusion PDF with asymmetric structure 
or a Gaussian mixture model. If one chooses to select one model rather than 
the other (i.e., choose an asymmetric and scalene PDF or multiple-tensor), 
then this decision is based more on the underlying assumptions of the model 
rather than on the evidence directly provided by the observed data. A large 
(possibly infinite) collection of Gaussian diffusion processes may be used to 
approximate an observed set of measurements to an arbitrary accuracy, but 
one has to consider the possibility that the information being fitted is noise 
instead of signal. We believe the rule of parsimony should be exercised at all 
times, and that summaries of orientational structure can be estimated and 
interpreted in (/-space rather than using (potentially) over-parameterized 
models. 

One potential extension to the methods proposed here would be to acquire 
multiple shells of a fixed radius in g-space instead of typical HARDI sampling 
[Wu and Alexander (2007); Khachaturian, Wisco and Tuch (2007)], that is, 
multiple-wavevector or hybrid imaging. In this case the test statistics are 



NONPARAMETRIC TESTING FOR HARDI 



35 



calculated for each shell, and then averaged across the different shells. The 
dominant orientation would be estimated by a weighted averaging of the es- 
timated dominant orientations for each shell, since its distribution depends 
on the SNR that is shell-dependent. Another possible acquisition method 
is diffusion spectrum imaging (DSI), corresponding to a Cartesian sampling 
of the characteristic function [Wedeen et al. (2005)]. It is more difficult to 
achieve the same directional resolution in DSI versus multiple-wavevector 
imaging, and so with realistic sampling times it may not be feasible to per- 
form the same analysis as outlined in this paper. However, other nonpara- 
metric summaries could be defined directly in g-space to characterize the 
spatial properties. 

One potential application of these g-space summaries would be to improve 
fiber-tracking algorithms, similar to the use of the Hessian of a local peak to 
improve probabilistic tractography models [Seunarine et al. (2007)]. These 
summaries would be used in addition to directions, to allow more careful 
tracking through forking and fanning structures (Figure 7), and distinguish 
local structure more consistently with crossing from such features using both 
the asymmetry and ellipsoidality measures. 
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